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Abstract 


1.  Background 

Improvement  in  the  detection  and  discrimination  of  UXO  using  the  EMI  regime  has  been  consis¬ 
tent  during  the  last  two  decades,  in  both  instrument  hardware  and  detection  software.  There  is 
nonetheless  still  a  disparity  between  the  often  idealized  environments  considered  (“golf  course” 
situation)  and  the  reality  in  the  field.  Working  on  diminishing  this  disparity,  various  SERDP 
programs  focus  on  developing  versatile  (next  generation)  instruments  as  well  as  robust  software 
for  data  processing  in  realistic  environments,  yielding  ever  improved  classification  performance. 

The  current  protocol  for  surveying  fields  is  typically  organized  in  two  stages.  First  the  sensors 
are  used  in  a  survey  mode,  or  dynamic  mode,  collecting  data  over  large  areas  and  flagging 
locations  which  may  be  contaminated  by  a  subsurface  UXO.  The  flagging  process  is  typically 
based  on  magnetic  field  amplitude,  based  on  the  rationale  that  high  amplitudes  may  correspond 
to  buried  UXO  and  should  therefore  be  further  investigated.  Hence,  at  the  second  stage,  the 
sensors  return  to  the  locations  of  interest  for  more  in  depth  data  collection  in  order  to  feed  the 
identification  and  classification  algorithms  with  high  quality  and  diverse  data.  Flagging  based  on 
field  amplitude,  however,  has  potential  drawbacks:  the  magnetic  field  in  the  EMI  regime  decays 
fast  and  a  deep  UXO  may  not  produce  a  strong  response  while  still  being  potentially  hazardous. 

2.  Objectives 

The  purpose  of  this  SEED  project  is  to  investigate  an  alternative  method  to  flag  areas  of  interest 
at  the  first  stage  of  data  collection  with  the  following  objectives: 

•  Work  with  dynamic  data  -  Dynamic  data  are  typically  more  noisy  and  less  diverse  than 
cued  interrogation  data  (collected  at  the  second  stage).  The  method  should  therefore  be 
tolerant  to  noise  and  able  to  provide  inversion  results  with  such  a  reduced  dataset. 

•  Provide  position  and  polarizability  estimates  -  We  wish  to  go  beyond  the  field  plot  and 
better  identify  the  targets,  in  order  to  make  more  informed  decisions  on  the  probability 
of  presence  of  targets  of  interest.  The  output  of  the  method  should  therefore  include 
target  position  in  the  (xy)  plane  for  subsequent  data  collection  if  necessary,  target  depth 
to  still  potentially  flag  deeper  UXO,  and  target  polarizabilities  to  have  a  first  identification 
procedure  and  not  miss  deep  UXO  or  small  items. 

•  Real-time  processing  -  The  method  should  not  slow  down  the  survey  speed  and  should 
therefore  work  in  real-time  while  still  satisfying  the  two  constraints  above.  For  purpose  of 
practicality,  we  set  our  real-time  limit  to  100  ms,  which  is  the  rate  at  which  dynamic  data 
are  currently  acquired  by  sensors  such  as  the  MPV-II  and  the  MetalMapper. 
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4.  Summary  of  results 


•  Applicability  to  small  items  -  The  method  should  be  validated  on  multiple  UXO  including 
some  challenging  ones  such  as  20-mm  projectiles. 

•  Validation  with  next  generation  sensors  -  Actual  data  from  dynamic  measurements  should 
be  used  as  the  ultimate  validation  of  the  method. 


3.  Summary  of  Science/Technology 

The  approach  proposed  by  this  work  is  to  implement  an  algorithm  based  on  Kalman  filters 
(KF)  and  Extended  Kalman  filters  (EKF)  for  the  sequential  processing  of  data  as  they  become 
available.  Kalman  filters  have  been  developed  as  a  stochastic  signal  processing  technique  using 
Bayes  rule  in  order  to  filter  a  signal  out  of  a  well  characterized  noise.  Kalman  filters,  however, 
only  provide  an  optimal  minimum  mean  square  error  estimator  if  some  constraints  are  satisfied. 
We  supposed  throughout  the  report  that  these  constraints  are  satisfied:  all  noise  sources  can  be 
drawn  from  a  Gaussian  distribution  as  well  as  prior  probability  density  functions  related  to  the 
unknown  parameters  to  be  estimated.  Kalman  filters  also  require  the  parameters  to  be  linear  with 
measurements:  within  the  dipole  approximation  often  used  to  process  EMI  data,  this  condition 
is  satisfied  by  the  dipole  moment  /77(f)  which  can  be  directly  related  to  the  polarizability  tensor 
/3(f).  The  position,  however,  appears  as  a  non-linear  parameter  which  we  estimate  using  the 
extended  Kalman  filter  and  a  linearization  of  the  equations.  Finally,  both  algorithms  (KF  and 
EKF)  are  iterated  to  sequentially  process  the  data  and  yield  converged  values  of  all  parameters 
to  be  estimated. 


4.  Summary  of  results 

The  main  achievement  of  this  SEED  project  is  the  implementation  and  validation  of  an  iterative 
KF-EKF  approach  satisfying  all  conditions  above:  dynamic  data  processing,  inversion  of  all  the 
parameters  of  the  dipole  model,  and  real-time  processing.  In  addition,  the  validation  has  been 
performed  on  various  targets,  including  a  small  20-mm  item  which  remains  a  challenging  target 
to  identify.  The  validation  of  the  KF-EKF  iterative  algorithm  has  been  performed  on  all  next 
generation  sensors  available  to  us:  TEMTADS,  MPV,  MPV-II,  and  MetalMapper.  For  practicality 
purposes,  we  have  concentrated  on  the  last  two  sensors  since  the  TEMTADS  is  usually  not  used 
for  dynamic  surveys  and  since  the  MPV-II  is  currently  replacing  the  MPV. 

It  is  worth  mentioning  that  the  purposes  of  dynamic  data  collections  with  the  MPV-II  and 
the  MetalMapper  are  very  different  in  our  investigation.  The  MPV-II,  meant  for  hand-held 
operation,  was  used  to  interrogated  a  limited  area  within  which  a  single  target  was  present, 
akin  to  being  waved  by  an  operator  in  a  detection  mode.  The  MetalMapper  instead  was  driven 
along  adjacent  60-m  long  lanes  to  flag  target  locations  across  a  large  geographical  area.  In  this 
latter  configuration,  multiple  targets  were  present  in  the  underground  and  sequentially  appeared 
and  disappeared  from  the  field  of  view  of  the  sensor  as  it  was  driven  across  the  area.  In  both 
cases,  the  sensors  operated  in  dynamic  mode,  meant  for  rapid  land  survey  at  the  expense  of  data 
quality.  Despite  the  noisy  data,  the  KF-EKF  algorithm  was  able  to  converge  to  proper  solutions 
as  verified  independently  either  using  our  Gauss-Newton  reference  algorithm,  using  ground  truth 
information  whenever  available,  or  comparing  with  independently  obtained  results.  In  addition, 
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5.  Conclusions  and  implications  for  future  research 


processing  times  of  each  new  data  collection  were  measured  to  be  less  than  100  ms  on  a  regular 
2  x  2.8  GHz  Quad-Core  Intel  Xeon  computer,  which  is  the  real-time  limit  for  our  application. 
The  inverted  parameters  include  position,  depth,  and  polarizabilities,  albeit  limited  to  the  short 
data  collection  time  range  inherent  to  a  dynamic  survey.  Yet,  even  if  this  dynamic  data  is  not 
sufficient  to  perform  classification,  it  is  often  sufficient  to  identify  the  subsurface  anomaly  as 
possible  UXO  or  not,  either  based  on  simple  signal  correlation  to  a  library  or  based  on  a  volume 
estimate.  The  method  can  therefore  be  used  in  land  survey  for  real-time  target  mapping,  flagging 
those  locations  to  which  the  instrument  needs  to  return  for  more  exhaustive,  cued  interrogation, 
data  collection. 

An  IEEE  Transactions  on  Geoscience  and  Remote  Sensing  article  has  been  submitted  as 
the  result  of  this  project:  “Real-time  processing  of  electromagnetic  induction  dynamic  data  for 
unexploded  ordnance  detection’ ,  by  Tomasz  M.  Grzegorczyk  and  Benjamin  E.  Barrowes. 


5.  Conclusions  and  implications  for  future  research 

As  part  of  this  work,  a  GUI  was  proposed  to  display  the  results  to  an  operator,  showing  in  real 
time  the  evolution  of  position,  depth,  and  polarizability  estimates  as  new  data  become  available 
during  the  measurement  process.  Combining  this  information,  an  indication  on  target  aspect 
ratio,  a  volume  estimate,  and  the  likelihood  of  being  a  target  of  interest  are  also  displayed  on 
the  GUI.  Such  an  integrated  tool  could  be  a  useful  feature  to  implement  on  existing  sensors  and 
taken  on-board  during  field  surveys.  Future  work  in  this  direction  would  include: 

1.  Integration  to  existing  on-board  softwares  with  sensors  of  interest 

2.  Extensive  validation  in  various  terrain  configurations 

3.  Optimization  of  feedback  information  (probability  of  target,  library  matching,  volume  esti¬ 
mation) 

4.  Extension  to  new  sensors  (in  particular  Pedemis  [Barrowes  12]) 

5.  Use  of  an  8  ms  dynamic  window  (instead  of  the  current  2.7  ms)  to  provide  considerably 
improved  identification  possibilities  since  the  KF-EKF  algorithm  achieves  more  than  just 
detection 
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The  detection  and  discrimination  of  unexploded  ordnance  (UXO)  remains  an  important  environ¬ 
mental  and  engineering  challenge  which  has  attracted  much  attention  for  the  last  two  decades. 
Accomplishments  reported  from  both  hardware  and  software  point  of  views  have  moved  the  re¬ 
search  focus  from  simple  laboratory  tests  to  much  more  complex  environments,  accompanied  by 
instruments  and  inversion  algorithms  adapted  to  ever  more  realistic  terrains.  This  trend  is  for  ex¬ 
ample  illustrated  by  the  recent  evolution  of  electromagnetic  induction  (EMI)  instruments  toward 
ever  improved  versatility  and  accuracy.  From  the  GEM  family  [Won  99]  to  TEMTADS  [Nel¬ 
son  07,  Nelson  01],  MPV  [Fernandez  11],  and  Metal  Mapper  [Prouty  09b],  the  instruments  have 
gained  properties  such  as  efficiency,  portability,  vectorial  excitation,  vectorial  reception,  etc.  The 
transmitter  and  receivers  are  often  mounted  on  the  same  sensor  and  are  therefore  co-located  [Nel¬ 
son  07,  Prouty  09b,  Fernandez  11],  although  recent  sensors  attempt  to  decouple  the  transmitter 
and  receiver  location  for  broader  coverage  and  more  diverse  information  content  [Barrowes  12], 
Cart-mounted  instruments  have  the  important  advantage  of  providing  consistent  measurements 
from  transmitter  to  transmitter,  and  can  measure  either  a  single  component  of  the  magnetic 
field  (e.g.  TEMTADS)  or  all  three  components  (e.g.  MetalMapper).  Their  disadvantage,  how¬ 
ever,  comes  from  the  fact  that  these  sensors  cannot  be  deployed  in  challenging  terrains  where 
vegetation  or  surface  state  are  non  trivial.  For  these  cases,  hand-held  instruments  such  as  the 
MPV  are  more  desirable  and  open  the  possibility  of  surveying  areas  of  difficult  access. 

Accompanying  hardware  improvements,  increasingly  more  sophisticated  algorithms  have 
been  developed.  Common  to  most  methods  is  the  inversion,  either  sequential  or  simultane¬ 
ous,  of  the  intrinsic  characteristics  of  the  targets  (polarizabilities)  and  of  their  position  7.  It  has 
been  recognized  that  the  inversion  of  7  is  more  challenging  and  often  requires  the  solution  of 
an  ill-conditioned  system.  Various  methods  have  been  proposed  to  address  this  problem,  from 
direct  optimization  based  on  analytical  methods  [Zhang  01],  model-based  approaches  [Miller  01], 
statistical  approaches  [Tantum  01],  expansions  into  spheroidal  and/or  ellipsoidal  modes  [Bar¬ 
rowes  04,  Chen  07,  Grzegorczyk  08],  to  more  direct  approaches  based  on  scalar  and  vector  poten¬ 
tials  [Shubitidze  08].  The  dipole  moment  rh,  on  the  other  hand,  is  easier  to  compute  once  the 
position  is  known,  and  can  be  obtained  either  at  every  time  channel  (or  every  frequency  if  such 
instruments  are  used)  or  following  a  parametrized  equation  [Pasion  07],  More  recent  works  have 
proposed  a  simultaneous  inversion  of  both  parameters,  yielding  methods  that  have  been  applied 
to  the  inversion  of  multiple  targets  present  in  the  field  of  view  of  the  sensor  [Song  09,  Grzegor¬ 
czyk  09,  Grzegorczyk  11],  akin  to  multi-UXO  configurations  or  to  UXO  in  the  presence  of  a  few 
but  large  clutter  items. 

The  fast  field  decay  in  the  EMI  regime  and  the  depths  at  which  UXO  are  typically  buried 
(usually  10  cm  or  more)  results  in  negligible  near-field  effects  at  the  receiver  level  so  that  the 
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secondary  field  is  well  approximated  by  that  of  a  dipole  [Jackson  99,  Bell  01,  Zhang  03]: 


H(r,  t )  = 


47r|  rl 


3 7(7-  m(t )) 


—  m 


47r|r| 


T=&  ~  1  )  ■ 


(1.1a) 

(1.1b) 


where  7  —  fRx  —  7'  is  the  distance  between  the  receiver  at  position  fRx  and  the  dipole  modeling 
the  UXO  at  position  7' .  The  dipole  moment  m(t)  is  related  to  the  intrinsic  polarizabilities  of  the 
target  P(t)  [Bell  01]  by  a  rotation  about  the  target  angles  R{6,(p )  weighted  by  the  primary  field 
at  location  7': 

m(t)  =  RT (6 ,(p)  j0(t)  R{9,(p)  Hpr(7Tx  -  7')  (1.2a) 

fpx(t)  0  0  \ 

where  /3(t)  =  I  0  Py(t)  0  (1.2b) 

V  0  0  pz(t)J 

In  order  to  simplify  the  notation,  we  omit  the  dependence  of  m(t)  on  the  transmitter  position  rTx 
and  write  /r?(fTx  —  7',  t )  =  m(t).  The  polarizabilities  P(t)  are  intrinsic  features  of  the  target,  inde¬ 
pendent  on  position  and  orientation,  and  can  therefore  be  used  for  classification  purposes.  Various 
techniques  have  been  proposed  to  estimate  the  parameters  (7',6,<p,  P(t)),  from  physics-based 
inversions  [Shubitidze  08],  deterministic  approaches  [Pasion  01b,  Pasion  07,  Grzegorczyk  11],  to 
statistical  approaches  [Collins  02,  Pulsipher  04,  Aliamiri  07],  All  these  methods  have  been  shown 
to  yield  accurate  parameter  estimation  in  various  configurations,  from  single  to  multi-target,  and 
can  therefore  be  used  to  feed  identification  methods  [Pasion  Ola.Tantum  01,  Shubitidze  10,Kap- 
pler  11,  Beran  11a,  Beran  lib,  Kass  12]  as  well  as  classification  methods  [Fernandez  10,  Hu  04] 
with  an  dig/no-dig  decision  with  a  quantifiable  degree  of  confidence. 

A  proper  identification  and  classification  require  high  quality  data  which  in  turn  require  a 
good  positioning  of  the  sensors  atop  the  UXOs.  The  data  collection  protocol  is  currently  based 
on  a  two-stage  process:  (1)  target  mapping,  whereby  the  sensors  are  used  in  a  detection  mode 
to  simply  flag  areas  of  interest  and  (2)  cued  interrogation  whereby  the  sensors  return  to  the 
areas  of  interest  for  more  thorough  data  collection.  At  the  first  step,  data  quality  is  sacrificed 
at  the  profit  of  speed,  with  the  purpose  of  surveying  large  areas  quickly.  Target  mapping  is 
currently  based  on  field  measurements,  identifying  regions  of  high  amplitudes  likely  to  correspond 
to  actual  targets.  However,  this  procedure  can  be  incompatible  with  field  situations  since  deep 
targets  could  produce  weak  measured  secondary  field  due  to  the  strong  decay  in  the  EMI  regime, 
yet  still  corresponding  to  hazardous  items.  Improved  methods  for  target  detection  are  therefore 
necessary,  and  proposing  one  such  method  is  the  purpose  of  the  present  work. 

The  desirable  features  of  a  detection  algorithm  are  to  provide  estimates  on  the  target  loca¬ 
tion  and  signature,  while  in  addition  being  fast  in  order  to  be  compatible  with  a  rapid  land  survey. 
A  method  based  on  correlation  to  a  sphere  signal  was  proposed  in  [George  11]  and  demonstrated 
to  successfully  estimate  the  location  and  depth  of  targets.  The  assumption,  however,  is  that 
the  three  principal  axes  of  a  target  respond  with  equal  amplitude  within  the  short  time  decay 
window  measured  during  a  dynamic  survey  (typically  2.7  ms).  Although  this  is  usually  the  case 
at  very  early  time,  the  assumption  quickly  breaks  down  at  later  time  channels.  In  this  work, 
we  return  to  the  tri-axial  dipole  approximation  of  Eq.  (1.1b)  and  propose  a  method  based  on 
an  iteration  between  Kalman  filters  (KF)  and  Extended  Kalman  filters  (EKF)  to  estimate  the 
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parameters  m(t)  -  hence  /3(t)  -  and  ( r',6,<p )  without  constraints  beyond  the  statistics  required 
by  KF.  In  the  following  chapters,  we  briefly  review  the  fundamentals  of  KF  and  EKF  and  apply 
them  to  both  static  and  dynamic  measured  data,  with  a  special  emphasis  on  the  MPV-II  and  the 
M  et  a  I M  a  p  pe  rse n  so rs . 
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2.  Kalman  filters:  background  and  validation 


The  notation  used  throughout  this  report  is  as  follows:  vectors  and  matrices  appear  with  a  bar  and 
double  bar,  respectively,  subscript  indicate  discrete  realizations  of  a  random  process  which  we  call 
instance  as  opposed  to  the  more  commonly  used  time,  in  order  not  to  introduce  possible  confusion 
with  the  actual  time  decay  of  the  magnetic  field  as  the  eddy  currents  diffuse  into  the  targets. 
Although  this  real  time  decay  could  be  used  within  the  Kalman  filter  framework  we  formulate 
subsequently  (because  of  the  assumption  of  time-invariant  UXO  position  and  orientation  angles), 
we  refrain  from  doing  so  simply  because  of  the  nature  of  the  sensors  we  use  which  provide  the 
entire  measured  time  decay  of  the  measured  magnetic  fields  all  at  once  and  not  sequentially. 

1.  Review  of  filters  implementation 

Kalman  filters  have  been  developed  as  a  stochastic  signal  processing  technique  using  Bayes'  rule 
in  order  to  filter  a  signal  out  of  a  well  characterized  noise.  Although  they  are  an  important 
generalization  over  Wiener  filters  in  that  they  can  accommodate  vector  signals  and  noises  which 
may  not  be  stationary  (which  are  properties  of  interest  to  us  in  the  realm  of  UXO  detection), 
Kalman  filters  only  provide  an  optimal  minimum  mean  square  error  estimator  if  some  constraints 
are  satisfied.  The  latter  are  stated  in  the  brief  review  included  hereafter,  understanding  that  this 
review  is  merely  meant  to  facilitates  the  comprehension  of  the  remainder  of  the  report  and  not  to 
provide  an  exhaustive  and  mathematically  complete  demonstration  of  the  underlying  derivations 
of  Kalman  and  Extended  Kalman  filters.  The  interested  reader  is  referred  instead  to  specific 
literature  on  this  topic  [Kalman  60,  Kay  93,  Gordon  02], 

We  start  by  considering  a  general  state  sequence  {xk,  k  e  N}  whose  evolution  is  written  as 


Xk  =  fk(xK- 1,  H-i), 


(2.1a) 


where  fk  is  the  governing  function.  In  our  case,  the  state  vector  xk  is  composed  of  the  unknowns 
to  be  determined,  namely  the  components  of  the  position  vector  or  of  the  polarizability  of  the 
UXO.  These  unknowns  are  related  to  the  measurements  by  the  function  hk\ 


Vk  =  hk(xk,  nk), 


(2.1b) 


the  measurements  being  the  magnetic  field  components  collected  by  EMI  sensors  across  trans¬ 
mitters  and  receivers.  In  Eqs.  (2.1),  {vk_1,  k  e  N}  and  {n,  k  e  N}  are  independent  and  identically 
distributed  process  noise  sequences  with  covariance  matrices  Qk-i  and  Rk,  respectively. 

The  data  fjk  thus  collected  constitute  sequential  information  identified  by  the  subscript  k. 
Our  purpose  is  to  estimate  the  state  xk  at  instance  k  knowing  the  measurements  up  to  instance 
k  and  the  state  xk_1  at  instance  ( k  —  1).  In  other  words,  we  need  to  estimate  the  conditional 
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probability  density  function  (pdf)  p(xk\Vi,k ),  where  {fjlk,  k  G  N}  represents  the  measurement 
vector.  This  pdf  can  typically  be  obtained  recursively  in  two  steps:  predict  and  update.  At  the 
prediction  stage,  the  pdf  p(xk-i\xi,k-i)  is  known  at  instance  (k  —  1)  and  the  prior  pdf  at  instance 
k  is  obtained  using  the  Chapman-Kolmogorov  equation  [Gordon  02] 


P(xk \Vi.k-i) 


dxk_x  p(xk\xk_l)  p(xk-i\vi,k-i)- 


(2.2) 


When  measurements  at  instance  k  become  available,  the  posterior  pdf  is  obtained  using  Bayes' 
rule: 


P(xk\Vi,k) 


P(Vz\xk)  pjXklTh.k-i) 

p(Vk\fji,k-i) 


(2.3) 


where  p(fjk\fii,k-i)  is  a  normalizing  constant.  Kalman  filters  are  based  on  the  assumption  that 
the  posterior  pdf  at  every  step  is  Gaussian  and  as  such,  parametrized  by  a  mean  and  a  covariance 
matrix.  This  assumption  holds  if  the  conditional  pdf  at  instance  (k  —  1)  is  Gaussian  (with 
covariance  matrix  Pk-i),  if  both  signal  and  measurement  noises  are  drawn  from  a  Gaussian 
distribution,  and  if  fk  and  hk  are  linear  functions  of  the  parameters.  Under  these  assumptions, 
Eqs.  (2.1)  are  rewritten  as: 


xk  =  Fkxk- 1  +  vk- 1 ,  (2.4a) 

fjk=  Hkxk  +  nk,  (2.4b) 

which  leads  to  the  Kalman  filter  (KF)  algorithm  [Kay  93,  Gordon  02]: 

xk\k—i  Fk  xk— i\k— 1>  (2.5a) 

Pk\k-i  —  Qk-i  +  FkPk-i\k-iF^ ,  (2.5b) 

Xk\k  =  Xk\k_ !  +  Rk(fik  -  Hkxk\k_i),  (2.5c) 

Pk\k  —  Pk\k-i  ~  KkHkPk\k-i,  (2.5d) 

where  Kk  —  Pk\k-1Hl(HkPk\k-iHl  +  Rk)^-  (2.5e) 


These  relations  are  directly  applicable  within  the  framework  of  the  dipole  approximation  of 
Eq.  (1.1b)  since  the  dipole  moment  m(t)  is  a  linear  function  of  the  magnetic  field.  Within  the 
same  approximation,  the  position  and  UXO  orientation  appear  as  non-linear  parameters  which 
therefore  cannot  be  estimated  using  Eqs.  (2.5)  directly.  Instead,  we  resort  to  the  Extended 
Kalman  filter  (EKF),  based  on  the  linearization  of  the  state  and  measurement  equations.  The 
extended  Kalman  algorithm  is  given  by  [Kay  93,  Ristic  04]: 


xk\k- 1  —  fk(xk-l\k-l)  • 

Pk\k-1  —  Qk-1  +  FkPk-i\k-iFj , 
xk\k  =  Xk\k- 1  +  KkiVk  -  hk(xk\k-l)), 
Pk\k  —  Pk\k-1  —  RkHkPk\k-li 


where  Fk  =  % 
ox 


x—xk-l\k-l 

~T 


H  =  — 
k  dx 


x=*k\k-i 

T  ,  d  \-l 


Kk=  Pk\k-iH'k(HkPk\k-iH'k  T  Rk) 


(2.6a) 

(2.6b) 

(2.6c) 

(2.6d) 

(2.6e) 

(2.6f) 
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In  our  implementation,  however,  the  usual  update  equation  (2.6d)  for  the  covariance  matrix 
Pk \k  is  replaced  by  its  Joseph-stabilized  version  [Bucy  68]  in  order  to  improve  numerical  stability 

Pk \k  =  M  ■  Pk\k- 1  ■ MT  +  Kk-Rk  ■  KTk  (2.7) 

M  =  f- Rk  ■  Rk. 


A  convergence  comparison  between  the  two  forms  is  illustrated  in  Figure  2.1  on  synthetic 
data,  where  the  stability  of  the  Joseph  form  over  the  standard  formulation  appears  clearly,  as 
demonstrated  by  the  much  faster  convergence  rate  of  the  filters. 


■1.5 - 1 - 1 - 1 - 1 - 1 

0  50  100  150  200  250 

Measurement  sequence 


(a)  Example  1  using  Eq.  (2.6d) 


(b)  Example  1  using  Eq.  (2.7) 


100  150 

Measurement  sequence 


250 


(c)  Example  2  using  Eq.  (2.6d) 


(d)  Example  2  using  Eq.  (2.7) 


Figure  2.1:  Evolution  of  the  estimated  positions  (x,y,z)  with  the  measurement  sequence  from 
synthetic  dynamic  MetalMapper  data.  Initial  guesses  for  examples  1  and  2  are  (—29,  0,  —12)  [cm] 
and  (—20,  —20,  —6)  [cm]  for  a  final  answer  at  (—34,  —2,  —49)  [cm]. 
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2.  Application  to  UXO  detection 

The  application  of  Kalman  filters  to  UXO  detection  is  somewhat  unusual  because  targets  are 
stationary,  which  is  to  say  that  xk  —  xk_i  for  the  estimation  of  {?' ,  9,<p).  The  prediction  stage  can 
therefore  be  omitted.  Subsequently,  the  two  algorithms,  KF  of  Eqs.  (2.5)  and  EKF  of  Eqs.  (2.6- 
2.7),  are  iteratively  applied  to  measured  data  in  order  to  sequentially  estimate  the  positions  and 
orientations  (r\  9,<p)  -  nonlinear  parameters  -  and  the  dipole  moment  m(t)  -  linear  parameters 
-  of  UXO.  The  algorithm  is  detailed  in  Table  2.2  and  reveals  the  iterative  procedure  between 
KF  and  EKF.  The  definitions  of  x,  f,  and  h  therefore  vary  with  the  parameters  that  need  to  be 
estimated.  The  dipole  moment  m  is  estimated  with: 


Xk  =  mk, 

(2.8a) 

Fk  =  l 

(2.8b) 

4= 

47T  T  3  V  \rk\2  ) 

(2.8c) 

Note  that  Fk  —  T  is  valid  only  at  a  given  transmitter  position,  which  is 

correct  within  step  (1) 

of  the  algorithm  of  Table  2.2.  The  estimation  of  (7',  9,<p)  supposes 
independent  with  the  following  definitions: 

these  parameters  time 

xk  =  rk,  fk{xk)  =  xk,  hk  =  Qr 

(2.8d) 

-  ,  -  r  dHk 

Xk  =  (pk,  fk{xk )  =  xk,  hk  = 

ocpk 

(28e) 

-a  ff-f  -  7  dRk 

xk  =  9k,  fk(xk)  =  xkl  hk  = 

ouk 

(2.8f) 

where  the  magnetic  field  Hk  is  obtained  from  the  dipole  model  of  Eq.  (1.1b).  In  practice,  the  three 
definitions  of  x  are  combined  into  a  single  vector  and  the  corresponding  hk  operator  is  obtained 
by  concatenating  the  Jacobian  matrices.  The  latter  require  the  derivatives  of  the  primary  field 
Hpr(r),  as  seen  in  Eq.  (1.2a),  which  can  be  computed  analytically  since  all  sensors  of  interest  to 
us  have  either  circular  or  square  primary  coils  [Stratton  41,  Jackson  99]. 

2. a  Polarizability  inversion 

In  order  to  validate  the  KF  algorithm,  Figure  2.2  shows  inversion  results  of  the  dipole  moment 
ifi  in  the  case  of  synthetic  data  with  and  without  noise,  when  the  position  of  the  dipole  is  known. 
Data  are  here  composed  of  the  z  components  of  the  magnetic  field  collected  over  a  grid  of  similar 
extent  to  the  TEMTADS  instrument  (square  grid  of  5  x  5  receivers  spaced  by  40  cm  in  both 
directions).  The  figures  show  the  evolution  of  the  estimated  parameters  (the  three  components 
of  m)  as  the  sequence  of  data  becomes  available  and  starting  from  a  poor  initial  guess.  It  is  seen 
that  the  filter  needs  about  seven  steps  to  lock  onto  the  right  values  in  a  noiseless  configuration.* 
In  a  very  noisy  configuration  and  assuming  that  noise  is  well  characterized,  Figure  2.2(b)  indicates 
that  the  filter  needs  more  measurements  but  is  still  able  to  rapidly  lock  onto  the  correct  values. 

*l\lote  that  this  results  is  merely  an  illustration  of  the  convergence  of  the  Kalman  filter.  In  practise,  the  linear 
part  can  be  directly  inverted  for  if  no  noise  corrupts  the  signal. 
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No  noise 


20%  max(Hz)  noise 


Figure  2.2:  Kalman  filter  prediction  of  a  single  dipole  for  synthetic  measurement  data:  (a) 
noiseless  case;  (b)  original  data  to  which  Gaussian  noise  has  been  added  of  zero  mean  and 
variance  equal  to  20%  of  the  maximum  value  of  the  magnetic  field.  The  true  values  of  the  dipole 
are  m  —  (2,  0.2,  0.2),  7  —  (10,  20,  —40)  cm. 


2.b  Position  inversion 

Unlike  the  KF  algorithm,  the  EKF  algorithm  has  no  optimality  properties  and  its  performance 
depends  on  the  accuracy  of  the  linearization  which  must  be  verified  in  our  situation.  Figure  2.3 
provides  such  verification:  results  have  been  obtained  using  the  Hz  field  from  a  single  dipole  at 
position  7  —  (10,20,-40)  cm  and  moment  m  —  (2,  0.2,  0.2)  computed  over  a  grid  of  5  x  5 
points  of  similar  extend  to  the  TEMTADS  sensor.  It  should  be  emphasized  that  this  result  is 
not  a  demonstration  but  a  mere  verification  of  the  applicability  of  the  extended  algorithm  to  our 
particular  case  (we  recall  that  the  extended  Kalman  filter  has  no  optimality  property).  Figure  2.4 
illustrates  the  propagation  of  the  pdf  through  the  operator,  based  on  TEMTADS  data  collected 
above  a  4.2-inch  mortar.  It  is  seen  that  starting  from  a  Gaussian  distribution,  the  final  pdf  is 
strongly  non-Gaussian,  as  expected.  However,  its  mean,  computed  at  7  —  (0,0,— 61.5)  cm,  is 
in  very  good  agreement  with  that  obtained  with  a  deterministic  method,  while  its  variance  is 
relatively  small. 


Results  obtained  by  the  EKF  algorithm  are  listed  in  the  last  column  of  Table  2.1  and  reveal 
a  few  centimeter  errors  compared  to  ground  truth.  Interestingly,  the  inversions  obtained  from 
TEMTADS  data  are  much  more  accurate  and  very  stable  as  well  as  robust  to  noise,  giving  good 
confidence  on  the  applicability  of  the  EKF  to  position  estimation. 
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2.  Application  to  UXO  detection 


Figure  2.3:  Kalman  filter  prediction  of  the  position  of  a  single  dipole  source  located  at  r  = 
(10,  20,  —40)  cm  and  of  dipole  moment  ni  —  (2,  0.2,  0.2). 


Sensor 

Configuration 

Ground  truth 

Inversion  based  on 
[Grzegorczyk  11] 

EKF  Inversion 

TEMTADS 

4.2”  mortar,  hor¬ 
izontal 

[0, 0,-60] 

[0, 1,-61] 

[3, 0,-63] 

Nose  piece,  verti- 

[0, 0,-27] 

[0,0, -30] 

[2, 0,-31] 

Ld  1 

Half  round,  hori¬ 
zontal 

[0, 0,-44] 

[0, 1,-48] 

[2, 1,-49] 

MPV 

81-mm 

[-23,22,-48] 

[-15,23,-48] 

[-10  18  -51] 

105-mm 

[-21,22,-60] 

[-13,23,-60] 

[-11  21  -62] 

BLU26 

[0,22,-34] 

[0,22,-34] 

[1  26  -38] 

57-mm 

[5,22,-50] 

[1,22,-50] 

[0  21  -55] 

60-mm 

[0,22,-46] 

[0,22,-48] 

[1  25  -50] 

MetalMapper 

Sphere 

1-st  quadrant 

[  28,  29,-36] 

[  25,  32,-44] 

Sphere 

2-nd  quadrant 

[-28,  23,-35] 

[-24,  25,-38] 

Sphere 

3-rd  quadrant 

[-26,-29,-35] 

[-22,-35,-39] 

Sphere 

4-th  quadrant 

[  29,-25,-36] 

[  38,-22,-39] 

Sphere 

centered 

[  0,  -1,-34] 

[  5,  -7,-38] 

Table  2.1:  Reference  configurations  for  various  sensors  and  inverted  positions  using  the  Gauss- 
Newton  method  of  [Grzegorczyk  11]  and  the  Extended  Kalman  Filter  (EKF)  where  polarizabilities 
are  known.  All  dimensions  are  given  in  centimeters. 
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2.  Application  to  UXO  detection 
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(a)  Initial  distribution. 
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(b)  Final  distribution. 


Figure  2.4:  Posterior  probability  density  function  (right)  obtained  from  the  initial  Gaussian  distri¬ 
bution  (left)  after  applying  the  extended  Kalman  filter.  The  final  pdf  is  non-Gaussian  as  expected. 
The  initial  and  final  means  and  variances  are:  r™an  =  (20,  0,  —30)  cm  and  omi  —  (0.1,  0.1,  0.1); 
Cet  =  (0,  0,  -61.5)  cm  and  oend  =  (0.02,  0.3,  0.01). 


2.c  Simultaneous  inversion  of  polarizabilities  and  position 

The  independent  inversions  of  position  and  polarizabilities  can  be  combined  into  an  iterative  al¬ 
gorithm  whereby  7  and  m  are  sequentially  updated  until  a  certain  level  of  convergence  is  reached 
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3.  Iterative  KF-EKF  algorithm 


(a)  Initial  state:  linearly  decaying  polarizabilities 
and  position  error  of  (80,  80,  60)  cm  in  three 
directions. 


.  -1  .  .0  .  _1 

10  10  10 

Times  [ms] 

(b)  State  after  10  iterations.  Position  error:  (1.5, 
0,  0.5)  cm. 


Figure  2.5:  Inversion  results  of  position  and  polarizabilities  when  sequentially  processed  by  an 
iterative  algorithm.  The  state  is  initially  far  from  the  ground  truth  and  presents  a  reasonable 
accuracy  after  10  iterations. 

(typically  measured  by  the  mismatch  between  the  computed  and  the  measured  magnetic  fields). 
The  algorithm  is  seeded  with  random  means  and  large  variances  for  both  parameters,  reflecting 
the  fact  that  we  have  little  knowledge  of  the  final  solution.  Figure  2.5  shows  how  such  algorithm 
performs  with  mortar  TEMTADS  data:  the  initial  position  is  far  from  the  truth  and  the  polariz¬ 
abilities  are  initialized  to  a  simple  linear  decay.  Figure  2.5(b)  shows  the  state  of  the  parameters 
after  10  iterations:  the  position  error  is  reduced  to  a  few  centimeters  in  all  three  directions  and 
the  polarizabilities  are  well  aligned  with  their  expected  values.  The  fast  convergence  of  the  algo¬ 
rithm  suggests  the  possibility  of  integrating  Kalman  filters  into  an  efficient  stand-alone  detection 
software  with  good  noise  tolerance. 

3.  Iterative  KF-EKF  algorithm 

As  described  in  Table  2.2,  the  two  sets  of  unknowns,  m(t)  and  (r',6,<p),  are  estimated  iteratively 
and  sequentially: 

1.  Iteratively:  The  KF  and  EKF  algorithms  estimate  independent  variables  that  serve  as  input 
to  each  other  (KF  estimates  m(t)  then  / 3(t )  which  is  an  input  to  EKF,  EKF  estimates 
(r7,  9,<fi )  which  is  an  input  to  KF).  The  two  algorithms  are  iteratively  applied  to  eventually 
yield  converged  values  of  both  unknowns.  The  sequence  used  to  iterate  between  KF  and 
EKF  is  the  sequence  of  transmitter  positions  throughout  the  data  acquisition  process:  as 
the  MPV-II  is  waved  atop  a  target  (see  Section  3)  or  as  the  MetalMapper  is  driven  along 
a  lane  (see  Section  4).  In  Table  2.2,  the  iterative  scheme  corresponds  to  sweeping  the 
sensor  position  rTx  and  computing  steps  (1)  and  (2). 

2.  Sequentially  (index  k  in  Eqs.  (2.5)  and  Eqs.  (2.6)):  Each  individual  algorithm  updates  the 
value  of  the  unknowns  by  sequentially  processing  the  data  as  they  become  available. 
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3.  Iterative  KF-EKF  algorithm 


—  Initial  guess  for  r',9,<p,  j 3(t ) 

—  Sweep  transmitter  positions  fTx 

(1)  KF:  estimate  /3(t) 

•  Compute  primary  field  at  r' 

•  Compute  m(t)  using  Eq.  (1.2a) 

•  Sweep  receiver  positions 

Update  m(t )  using  Eqs.  (2.5) 

End 

•  Compute  (3(t)  using  Eq.  (2.9) 

(2)  EKF:  estimate  (r',0,0) 

•  Sweep  receiver  positions  fRx 

Sweep  time  channels 

Update  (r',0,0)  using  Eqs.  (2.6) 

End 

End 

-  End 

Table  2.2:  KF-EKF  iterative  algorithm  for  the  estimation  of  UXO  position  r' ,  orientation  (6,<p), 
and  time  dependent  diagonal  polarizability  tensor  /3(t). 


(a)  KF:  the  dipole  moment  m(t )  is  estimated  by  sweeping  all  receivers  sequentially  (at  a 
given  fTx  position).  The  estimation  is  performed  at  each  time  channel,  from  which 
the  time-dependent  polarizabilities  are  obtained  by  rewriting  Eq.  (1.2a)  as: 

R(d,4>)-  rh(t)  —  /3(t)  •  R(6,<l>)-  Hpr.  (2.9) 

(b)  EKF:  since  (f',0,0)  are  supposed  time  independent,  the  sequential  information  to  up¬ 
date  these  parameters  is  constituted  by  the  sequence  of  receivers  and  the  sequence  of 
time  channels.  Hence,  a  larger  dataset  is  available  to  estimate  the  nonlinear  parameter, 
which  is  a  desirable  feature. 

The  sequential  processing  in  the  algorithm  of  Table  2.2  corresponds  to  the  application  of 
Eqs.  (2.5)  at  every  time  channel  by  sweeping  the  sensor  positions  rRx,  and  Eqs.  (2.6)  by 
sweeping  both  time  channels  and  sensor  positions  FRx. 

In  the  subsequent  two  chapters,  the  algorithm  of  Table  2.2  is  validated  on  data  collected 
by  the  MPV-II  and  the  MetalMapper  sensors.  In  the  first  case,  we  consider  both  static  and 
dynamic  data  collected  on  a  grid  over  a  well-controlled  target.  In  the  second  case,  dynamic 
data  are  processed  from  large  geographical  area  collections  at  Camps  San  Luis  Obispo  (CA)  and 
Butner  (NC). 
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3.  Validation  with  the  MPV-II  sensor 


1.  The  sensor 

The  MPV-II  (Man  Portable  Vectorial  sensor  II)  is  a  state  of  the  art  EMI  sensor  for  the  detection 
and  classification  of  UXO,  bringing  important  flexibility  and  accuracy  improvements  compared 
to  the  first  generation  of  sensors  such  as  the  EM-63  [Ware  03]  for  example.  Some  of  the  key 
features  of  the  MPV-II  include  its  updated  electronics  and  its  vectorial  receivers  that  capture 
the  three  components  of  the  secondary  magnetic  fields.  In  addition,  the  MPV-II  is  compact 
and  lightweight  as  shown  in  Figure  3.1,  meant  to  be  operated  in  a  hand-held  mode  by  a  single 
operator.  The  unique  horizontal  transmit  coil,  of  diameter  49.68  cm,  illuminates  a  small  region 
and  needs  to  be  waved  atop  an  area,  either  randomly  or  over  a  pre-defined  grid,  to  collect  a 
sufficiently  large  amount  of  data  for  UXO  detection  and  classification. 

The  MPV-II  can  be  operated  in  two  modes,  referred  to  as  static  (or  cued  interrogation) 
and  dynamic.  In  this  report,  we  focus  on  the  dynamic  mode,  meant  for  rapid  land  survey:  the 
data  collection  window  is  limited  to  2.7  ms  instead  of  the  usual  25  ms  in  cued  interrogation 
mode.  Target  classification  is  therefore  practically  excuded  and  the  sensor  is  mostly  used  for 
detection  and  target  volume  estimation.  In  addition,  we  consider  the  MPV-II  to  be  waved  atop 
a  limit  spatial  region  of  about  0.25  m2,  thus  presumably  looking  for  a  single  target  only.  This 
is  in  contrast  with  the  large  survey  operation  mode  of  the  MetalMapper  discussed  subsequently 
where  targets  would  come  and  go  as  the  sensor  is  driven  along  a  60  m  long  lane. 

In  terms  of  modelling,  the  MPV-II  is  very  similar  to  the  first  generation  MPV:  circular 
transmit  loop  and  cubic  vectorial  receivers.  The  main  differences  in  configuration  are  the  physical 
locations  of  the  receivers  and  the  number  and  positions  of  the  loops  [Fernandez  11],  which  are 
all  differences  straightforward  to  account  for  in  modeling.  Consequently,  we  do  not  repeat  here 
the  formulae  implemented  to  compute  the  primary  field  and  its  derivative,  but  instead  refer  the 
reader  to  [Grzegorczyk  11], 


2.  Dynamic  mode 

There  are  several  differences  between  a  dynamic  and  a  static  mode,  one  of  them  being  the 
reduced  number  of  averaging  cycles  used  during  dynamic  data  collection  compared  to  a  static 
data  collection.  As  a  result,  the  signal  to  noise  ratio  (SNR)  is  typically  degraded  in  dynamic  data, 
resulting  in  non-smooth  field  decay  curves.  This  phenomenon  is  illustrated  in  Figure  3.2  across 
all  five  receivers  and  all  three  field  components.  The  figure  shows  the  field  values  measured  at 
various  cycles,  along  with  their  averaged  values  (black  overlapping  lines).  It  is  seen  for  example 
that  receiver  1  essentially  measures  only  noise  and  that  in  general  the  Hx  components  is  more 
noisy  than  the  Hy  or  Hz  component  (for  this  specific  example).  The  data  input  to  our  KF-EKF 
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2.  Dynamic  mode 


(a)  Overview  of  the  sensor.  (b)  Zoom  on  sensor  head. 


Figure  3.1:  Photographs  of  the  MPV-II  sensor,  comprising  a  boom  for  hand-held  maneuverability 
at  the  end  of  which  is  located  the  EMI  head.  The  latter  is  composed  of  a  circular  transmitter 
(49.68  cm  diameter)  and  five  equi-spaced  cubic  receivers  (8  cm  side). 


algorithm  consisted  of  the  averaged  field  values  as  shown  by  the  black  curves:  typically  smoother 
than  the  curves  at  each  individual  cycle,  but  more  noisy  than  a  measurement  from  a  static  data 
collection. 
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2.  Dynamic  mode 
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Figure  3.2:  Field  values  and  their  average  (solid  black  line)  from  a  series  of  dynamic  measurements 
centered  over  a  20  mm  target. 
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3.  Results 


3.  Results 

Reference  measurement  data  were  available  for  three  targets:  a  60-mm,  a  57-mm,  and  a  20-mm, 
all  shown  in  the  insets  of  Figure  3.3.  These  datasets  were  acquired  in  static  mode  and  are 
subsequently  used  to  simulate  dynamic  data  with  two  caveats.  First,  data  include  the  entire 
25  ms  window  of  magnetic  field  decay:  we  manually  truncate  the  dataset  to  the  first  2.7  ms. 
Second,  more  measurement  cycles  are  used  to  average  the  magnetic  fields  in  a  static  mode 
(typically  57  cycles  of  100  ms)  than  in  a  dynamic  mode  (as  few  as  10  to  25  cycles  are  used), 
resulting  in  a  higher  signal-to-noise  ratio  (SNR)  at  the  expense  of  a  slower  data  acquisition 
rate.  In  order  to  illustrate  the  effects  of  a  lower  SNR  on  inversion  results,  we  later  analyze  a 
dataset  collected  in  a  real  dynamic  mode.  For  this  first  step,  however,  we  resort  to  the  truncated 
static  data  collected  by  all  vectorial  receivers,  at  sensor  positions  spanning  a  nine-point  square 
grid  with  20  cm  point-to-point  separation.  In  all  cases  the  targets  were  located  under  the 
central  grid  point  and  measurements  were  started  at  the  top-left  grid  point.  The  polarizability 
curves  of  Figure  3.3  were  obtained  by  iterating  the  KF  and  EKF  algorithms  (colored  circles 
in  the  figures)  as  shown  in  Table  2.2  where  fTx  sweeps  the  nine-point  grid,  fRx  sweeps  the 
five  receivers,  and  the  time  channels  span  the  25  ms  of  data  collection.  The  reference  curves, 
in  gray  in  the  figures,  correspond  to  polarizabilities  obtained  with  our  Gauss-Newton  reference 
algorithm  [Grzegorczyk  11]  and  show  the  good  agreement  between  the  two  methods.  The  2.7  ms 
limit  is  indicated  for  information  purpose,  illustrating  the  fact  that  most  of  the  decay  information 
is  included  at  later  times,  unavailable  in  dynamic  datasets.  Consequently,  dynamic  data  are 
mostly  used  for  detection  and  volume  estimation  purposes,  whereas  classification  must  be  done 
by  returning  to  the  target  location  and  operating  in  cued  interrogation  mode  in  order  to  invert 
and  analyze  the  entire  time  decay  of  the  polarizabilities.  In  the  next  examples,  the  datasets  of 
Figure  3.3  are  truncated  to  the  first  2.7  ms  and  are  used  as  input  to  the  KF-EKF  algorithm  in 
order  to  simulate  a  dynamic  operation. 

We  first  study  the  position  convergence  capability  of  the  algorithm  as  function  of  initial  guess 
(7',  9,(p,  1 3(t ))  of  Table  2.2.  In  all  cases,  the  angles  are  started  at  zero,  the  polarizabilities  are  set 
to  be  equal  along  the  three  axes  with  a  1/t  slope  (where  t  denotes  time)  starting  from  an  arbitrary 
value  of  one  at  the  first  time  channel,  and  7'  randomly  samples  a  surface  area  between  ±1  m  in 
the  (xy)  plane,  and  set  to  a  constant  —40  cm  in  depth,  with  targets  ground  truth  position  at 
about  (20,  —20,  —30)  cm  as  indicated  by  the  green  circles  in  Figure  3.4.  The  algorithm  is  iterated 
across  the  nine-point  grid  for  the  transmitter  positions,  across  the  five  receivers  each  measuring 
three  components  of  the  magnetic  field,  and  across  30  time  gates  between  about  150_/us  and 
2.7  ms.  The  noise  covariance  matrices  are  set  to  Qk-i  —  Rk  —  T  x  10“20  while  Pk  —  T  x  10-1 
since  we  assume  no  a  priori  knowledge  of  true  target  position  and  therefore  possibly  a  bad  initial 
guess.  The  final  position  estimates,  indicated  by  the  blue  crosses,  all  end  up  within  5  cm  of  the 
ground  truth  (red  square)  and  therefore  indicate  a  good  convergence  of  the  algorithm  with  the 
available  quasi-dynamic  dataset.  By  rerunning  the  entire  KF-EKF  algorithm  with  initial  guesses 
between  ±50  cm,  the  converged  values  end  up  within  2.5  cm  of  the  ground  truth,  as  illustrated  in 
Figure  3.5  (with  similar  convergence  bounds  on  the  depth,  not  shown  here).  These  accurate  and 
consistent  position  estimates  reveal  that  the  2.7  ms  time  window  of  data  collection  is  sufficient 
to  properly  estimate  the  location  of  all  the  three  targets  under  investigation,  including  that  of 
the  small  20-mm  item. 
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3.  Results 
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Figure  3.3:  Polarizability  of  three  UXO  (see  insets)  from  static  data  measured  by  the  MPV-II. 
The  2.7  ms  limit  of  dynamic  data  collection  is  indicated.  The  reference  gray  signatures  are  based 
on  a  Gauss-Newton  inversion  [Grzegorczyk  11], 


Next  we  consider  an  actual  dynamic  dataset  collected  atop  the  small  20-mm  UXO.  Com¬ 
pared  to  the  truncated  static  dataset  of  the  previous  example,  the  signal-to-noise  ratio  of  dynamic 
measurements  is  deteriorated  owing  to  the  fewer  cycles  used  to  average  the  data.  A  comparison 
of  measured  field  values  in  static  and  dynamic  mode  is  presented  in  Figure  3.6.  The  configuration 
is  that  of  an  offset  target  compared  to  the  MPV-II  as  shown  in  Figure  3.6(a),  with  the  target 
along  the  y  direction.  Field  plots  at  the  closest  three  receivers  illustrate  the  smoother  nature 
of  the  static  dataset  in  all  three  components  of  the  measured  field.  In  particular,  the  noise  in 
the  Hz  component  (usually  dominant)  has  an  important  influence  on  the  inversion  results  and 
algorithm  convergence. 

The  KF-EKF  algorithm  is  run  in  a  similar  way  as  in  the  previous  example,  in  particular  with 
the  same  initial  values  of  the  covariance  matrices  Qk-i,  Rk.  and  Pk.  Since  the  algorithm  is  se- 
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(c)  60  mm  (d)  Zoom  on  (b) 


Figure  3.4:  Evolution  of  position  estimates  in  the  (xy)  plane  as  the  filter  is  updated  for  various 
transmitter  positions.  The  random  initial  guesses,  between  —1  m  and  +1  m  in  the  (xy)  plane, 
are  indicated  by  the  green  circles  and  the  final  estimated  positions  are  indicated  by  the  blue 
crosses.  The  ground  truths  are  indicated  by  red  squares.  The  gray  paths  illustrate  the  sequential 
estimates  during  the  iterative  process. 
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(c)  Zoom  on  (c) 

Figure  3.5:  Same  as  Figure  3.4  with  initial  guesses  between  —0.5  m  and  +0.5  m. 


quential,  however,  it  is  possible  that  one  noisy  measurement  at  a  given  instance  k  destabilizes  the 
solution  and  yields  a  divergence  (such  phenomenon  is  less  likely  to  happen  with  a  non-sequential 
algorithm  where  all  the  data  are  processed  simultaneously,  thus  minimizing  the  influence  of  one 
bad  measurement).  In  order  to  prevent  such  divergence  from  happening,  we  implement  a  memory 
scheme  whereby  the  newly  estimated  parameters  (r'k,  0k,  (pk,  mk(t ))  are  retained  only  if  they  yield 
a  lower  field  mismatch  than  the  previous  estimates  (rk_1, 6k-1,  4>k- 1,  mk-1(t)).  If  not,  the  (k  —  1) 
parameters  are  propagated  to  instance  k  and  the  KF-EKF  algorithm  is  iterated  to  the  ( k  +  1) 
instance.  The  polarizabilities  resulting  from  this  updated/select  scheme  for  the  20-mm  UXO  are 
shown  in  Figure  3.7  inverted  at  a  position  (—3,  —3,  —40)  cm  along  the  y  direction.  This  repre¬ 
sents  a  difference  of  3  cm  and  2  cm  in  the  x  and  y  directions,  respectively,  and  a  similar  depth 
to  the  inverted  position  using  the  truncated  synthetic  data  (UXO  inverted  at  (0,  —1,  —40)  cm). 
The  main  polarizability  in  Figure  3.7  follows  the  expected  decay  and  amplitude  properly  but  the 
two  weaker  ones  are  considerably  more  noisy,  which  we  attribute  to  the  lower  SNR  as  discussed 
above.  Nonetheless,  this  example  illustrates  the  possibility  of  using  the  MPV-II  sensor  for  the 
dynamic  detection  of  targets  as  small  as  the  20-mm. 
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(a)  Measurement  configuration 


Time  [ms] 

(b)  Receiver  3  (central) 


„„-7  H  field  at  receiver  4 
x  10 


(c)  Receiver  4 


(d)  Receiver  5 


Figure  3.6:  Comparison  between  field  components  collected  atop  an  off-centered  20-mm  UXO 
in  a  static  mode  (solid  lines)  truncated  at  2.7  ms  and  a  dynamic  mode  (dashed  lines). 


4.  An  alternative  formulation 


The  purpose  of  the  algorithm  presented  so  far  is  to  invert  for  r  and  P  in  Eq.  (1.1)-(1.2),  which  is 
an  approach  adopted  by  the  vast  majority  of  researchers  working  with  the  dipole  model.  Within 
the  KF-EKF  algorithm  developed  here,  P  are  inverted  using  a  KF  approach  (linear  part  of  the 
model)  while  r  is  inverted  using  an  EKF  approach  (non-linear  part).  In  general,  the  EKF  algorithm 
does  not  guarantee  convergence  and  its  efficacy  usually  depends  on  how  strongly  nonlinear  the 
problem  is.  In  order  to  ensure  a  better  convergence  of  the  EKF  algorithm,  it  is  possible  to  reduce 
the  nonlinearity  of  the  problem  by  inverting  for  modified  variables. 
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Figure  3.7:  Polarizabilities  of  a  20-mm  UXO  obtained  from  dynamic  MPV-II  measurements.  The 
reference  gray  signatures  are  based  on  a  Gauss-Newton  inversion  [Grzegorczyk  11], 


The  formulation  starts  by  writing  r  as 

f  =r  r  =  rTx+rRx-r'  (3.  la) 

where  r—  xsin  0  cos0  +  y  sin  0  sin  0  +  zcos0  (3.1b) 

where  r  is  the  same  distance  vector  as  in  Eqs.  (1.1). 

In  general,  the  angles  (0,0 )  are  not  good  parameters  to  invert  for  since  they  depend  on  both 
the  transmitter  rTx  and  receiver  7Rx  positions  (in  addition  to  depending  on  the  position  of  the 
UXO  r').  Fience,  an  inversion  algorithm  that  collectively  uses  all  data  (from  all  Tx  and  all  Rx) 
gains  little  by  adopting  (0,0)  as  unknowns  since  the  latter  vary  at  each  data  point.  The  situation 
is  very  different  for  a  sequential  estimation  algorithm  such  as  the  KF-EKF  algorithm:  (0,0)  can 
be  adjusted  by  the  known  quantities  rTx  and  7Rx  at  every  iteration  of  the  algorithm,  leaving  only 
local  angles  to  be  inverted  for  in  order  to  estimate  ?' .  The  non-linearity  of  the  magnetic  field 
with  respect  the  unknowns  (0,0)  is  much  reduced  compared  to  the  non-linearity  with  respect  to 
r. 

The  angles  (0,0)  provide  information  on  the  direction  at  which  the  UXO  is  located,  and  the 
remaining  parameter  1/r3  provides  the  range  information,  i.e.  how  far  in  the  given  direction  is 
the  UXO.  Flere  again,  in  order  to  limit  the  nonlinearity  of  the  magnetic  field  with  the  unknown, 
we  solve  for  1/r3  instead  of  solving  for  r  itself,  contributing  to  more  stability  of  the  EKF  algo¬ 
rithm.  We  also  note  that  the  dyad  r?  as  well  as  its  derivative  with  respect  to  0  and  0  are  all 
straightforward  to  compute. 

Consequently,  we  can  reformulate  the  algorithm  to  solve  for  (0,0,  1/r3)  and  process  the 
data  using  the  EKF  approach.  The  polarizabilities,  remaining  linear,  are  still  obtained  from  the 
KF  approach. 

An  example  of  inversion  is  presented  in  Figure  3.8,  comparing  the  predictions  between  the 
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Figure  3.8:  Polarizabilities  obtained  by  an  iterative  KF-EKF  sequential  algorithm  from  MPV-II 
data  collected  over  a  57-mm  UXO.  The  reference  gray  signatures  are  based  on  a  Gauss-Newton 
inversion  [Grzegorczyk  11],  The  inverted  positions  are  (4, —5, —42)  cm  and  (0,  — 3,  — 41)  cm 
obtained  with  the  KF-EKF  and  GN  algorithms,  respectively,  whereas  the  ground  truth  is  listed 
as  (0,  0,  —35)  cm  to  the  tip  of  the  vertical-nose-up  57-mm. 


KF-EKF  algorithm  and  our  Gauss-Newton  reference,  from  data  collected  over  a  57-mm  UXO. 
The  robustness  of  the  results  to  initial  guess  and  noise  levels  (not  shown  here)  confirms  the 
improved  stability  of  the  modified  algorithm  and  proposes  it  as  a  possible  alternative  formulation 
to  investigate  in  the  future. 
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1.  The  sensor 

The  MetalMapper  is  also  considered  a  next  generation  sensor  like  the  MPV-II.  Among  its  unique 
features  are  three  large  perpendicular  transmit  coils  (1  m  x  1  m)  with  electronically  switched 
TEM  transmitter  loop  drivers  that  provide  illumination  diversity  of  the  targets  from  a  single 
sensor  position,  and  seven  fully  vectorial  receivers  similar  to  those  of  the  MPV-II.  Additional 
hardware  and  operational  information  can  be  found  in  [Prouty  11],  A  photograph  and  a  schematic 
representation  of  the  receiver  locations  are  shown  in  Figure  4.1:  the  two  perpendicular  vertical 
transmitter  coils  are  clearly  visible  whereas  the  horizontal  transmitter  is  enclosed  in  a  wooden 
casing  along  with  the  receivers  (hence  invisible).  Like  the  MPV-II,  the  MetalMapper  can  be 
operated  in  either  a  static  -  or  cued  interrogation  -  mode  whereby  full  data  sets  are  collected 
atop  previously  flagged  locations  and  used  for  classification  purposes,  or  in  a  dynamic  -  or  target 
mapping  -  mode,  whereby  data  are  acquired  at  a  fast  rate  while  the  sensor  is  moving.  In 
2009  and  2010  the  sensor  participated  in  live  site  demonstrations  at  the  former  Camp  San  Luis 
Obispo  California  (SLO)  and  the  former  Camp  Butner,  North  Carolina,  in  order  to  validate  both 
operation  mode  in  realistic  field  situations.  In  a  cued  interrogation  setting,  the  three  transmit 
coils  fire  sequentially  to  excite  independent  eddy  currents  in  the  subsurface  targets,  and  the 
secondary  magnetic  field  components  (Hx,  Hy,  Hz )  are  collected  by  the  seven  receivers  over  a 
decay  period  of  about  25  ms.  The  dataset  is  therefore  composed  of  63  time  decay  magnetic 
field  curves  (three  transmitters,  seven  receivers,  three  components)  which  need  to  be  processed 
and  converted  into  position  and  time-dependent  intrinsic  polarizabilities  of  targets.  Various  field 
tests  have  shown  that  the  dataset  collected  by  the  MetalMapper  is  compatible  with  real  field 
detection  and  classification  of  various  UXO  [Prouty  09b,  Prouty  09a], 

In  this  report,  however,  we  are  interested  in  the  dynamic  mode,  meant  for  rapid  land  surveys 
and  for  the  creation  of  target  maps  to  which  operators  can  return  if  necessary.  The  main  feature 
of  a  dynamic  operation  is  therefore  speed.  At  the  data  acquisition  stage,  it  is  achieved  by  (1) 
using  only  the  horizontal  transmit  coil  (hence  one  transmitter  instead  of  three),  and  (2)  limiting 
the  recording  time  to  the  first  2.7  ms  of  field  decay  (instead  of  the  traditional  25  ms)  and  using 
less  averaging  windows,  thus  resulting  in  lower  signal  to  noise  ratios  than  in  cued  interrogation 
mode. 

The  protocol  at  the  two  demonstration  sites  mentioned  above  was  to  drive  the  MetalMapper 
along  lanes  separated  by  about  0.75  m  and  acquire  data  points  at  100  ms  intervals,  thus  creating 
magnetic  field  maps  of  the  area  to  be  surveyed.  Tracking  of  the  sensor  positions  was  performed 
through  an  on-board  RTK  GPS  (real-time  kinematic  global  positioning  system),  providing  a 
accuracy  of  about  2  cm  in  latitude  and  longitude  directions.  The  conversion  into  Northing  and 
Easting  coordinates  necessary  for  our  algorithm  was  performed  following  [Plager  89,  Smith  90]. 
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The  attitude  of  the  sensor  platform  was  obtained  through  a  two-axis  accelerometer  to  measure 
pitch  and  roll  angles,  and  a  three-axis  magnetometer  to  measure  heading  angle  with  reference  to 
the  true  magnetic  north  [Prouty  11],  This  information  is  necessary  in  order  to  properly  rotate  the 
magnetic  field  components  so  as  to  produce  an  estimated  position  in  a  global  reference  frame 
(arbitrarily  chosen). 

Inherent  to  purpose  of  land  survey  where  target  location  is  unknown,  every  dynamic  shot  is 
used  semi-independently,  i.e.  data  corresponding  to  successive  sensor  positions  are  not  combined 
into  a  single  dataset  but  are  used  independently.  However,  the  estimated  parameters  from  one 
sensor  position  are  used  as  input  at  the  next  position.  If  the  dominant  signal  between  two 
consecutive  sensor  positions  correspond  to  the  same  target,  the  initial  guess  is  good  and  helps 
the  filter  converge  more  rapidly.  If  not  (for  example  when  a  target  is  being  lost  because  too  far 
from  the  sensor  already,  or  when  a  second  target  produces  a  stronger  signal),  the  filter  needs  to 
adjust  the  estimated  parameters  from  the  inaccurate  initial  guess  to  the  new  solution.  The  data 
processing  follows  that  of  Table  2.2  exactly  where  the  outermost  loop  (sweeping  rTx)  corresponds 
to  the  sequence  of  transmitter  positions  along  each  lane. 


2.  Results  from  KF-EKF 

We  first  consider  the  dynamic  data  collected  at  Camp  Butner.  The  area  surveyed  is  approximately 
five  meters  by  60  meters,  sampled  into  eight  lanes  separated  by  about  75  cm  and  indicated  as 
thin  black  lines  in  Figure  4.2.  The  background  in  the  figure  corresponds  to  the  interpolated 
z-component  of  the  magnetic  field  at  the  central  receiver,  plotted  on  a  logarithmic  scale.  It  can 
be  seen  that  the  map  contains  many  peaks,  typical  of  a  soil  highly  contaminated  with  clutter 
items.  A  series  of  seeded  items  are  placed  along  a  line  at  an  Easting  relative  coordinate  of  about 
2.5  m.  Referring  to  the  figure,  the  items  are  (from  bottom  to  top)  a  shotput,  a  37-mm  projectile, 
three  small  ISO,  a  shotput,  and  a  pit,  all  identified  by  red  circles. 

The  selection  of  targets  from  inversion  results  is  currently  performed  on  a  per-lane  basis 
and  in  two  stages.  First  the  KF-EKF  algorithm  processes  every  data  point  semi-independently 
to  produce  an  estimated  target  position  and  signature,  regardless  of  the  fact  whether  there  is  a 
target  or  not.  Second,  the  estimated  positions  are  post-filtered  to  identify  those  detection  results 
that  correspond  to  actual  targets  based  on  result  consistency  in  Northing,  Easting,  and  depth 
inversions.  This  process  is  best  illustrated  with  an  example,  for  which  we  refer  to  Figure  4.3.  The 
figure  shows  the  sequence  of  positions  estimated  by  the  KF-EKF  algorithm  when  the  MetalMap- 
per  is  driven  along  lane  5,  at  relative  Northing  coordinates  between  30  m  and  32.5  m,  where  x 
refers  to  Easting,  y  to  Northing,  and  z  to  depth.  In  the  absence  of  target,  the  inverted  (x,y,z) 
coordinate  appear  random  from  one  sensor  position  to  the  next,  not  exhibiting  any  consistency 
or  relation  to  one  another,  allowing  the  post-processing  algorithm  to  filter  these  positions  out. 
In  the  presence  of  a  target,  however,  the  inverted  (x,y,z)  coordinates  exhibit  a  pattern  easily 
identifiable.  In  the  example  of  Figure  4.3,  the  strongest  signal  is  produced  by  the  last  shotput 
located  at  a  relative  Northing  coordinate  of  about  31.5  m.  The  inverted  y  coordinate  corresponds 
to  the  target  position  along  the  driving  direction  of  the  sensor:  the  target  is  first  in  front  of  the 
sensor  (positive  y),  then  at  the  same  Northing  level  as  the  sensor  (y  close  to  zero),  then  behind 
the  sensor  (negative  y).  At  30.5  m,  the  target  is  still  too  far  away  from  the  sensor  to  produce 
a  stable  inversion  and  the  coordinates  are  seen  to  be  quickly  varying,  indicating  no  convergence. 
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(a)  MetalMapper  sensor. 


(b)  Receivers  and  z  oriented  coil. 

Figure  4.1:  Geometry  of  the  MetalMapper  used  in  the  algorithm.  The  primary  field  coil  di¬ 
mensions  are  1  m  x  1  m  and  produce  a  magnetic  field  computed  from  the  addition  of  four 
wires  supporting  the  same  constant  current.  For  computational  efficiency,  the  received  field  is 
computed  at  the  center  of  each  receiver  only,  i.e.  the  finite  size  of  the  receiving  cubes  is  not 
accounted  for. 


Between  about  31  m  and  32  m,  the  inverted  y  position  is  seen  to  follow  the  dashed  black  line 
which  corresponds  to  the  local  y  position  of  the  sensor  (with  arbitrary  center):  a  good  match 
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indicates  a  stationary  target,  as  expected.  Within  this  range,  the  inverted  Easting  and  depth 
present  variations  due  to  the  noisy  data,  yet  the  amplitudes  of  these  variations  are  smaller  than 
when  the  algorithm  does  not  converge,  indicating  a  consistent  Easting  position  and  a  constant 
depth,  again  corroborating  the  detection  of  a  stationary  target.  Consequently,  inversion  results 
are  either  discarded  (around  32.5  m  for  example,  we  conclude  that  no  target  is  detected)  or 
flagged  as  targets  of  interest  based  on  the  pattern  in  (x,y,  z)  evolution  with  the  sensor  motion. 

Applied  to  Camp  Butner  data,  this  process  produces  the  estimated  positions  marked  by 
white  crosses  in  Figure  4.2  which  are  seen  to  usually  overlap  well  with  the  magnetic  field  peaks. 
In  some  cases,  however,  the  overlap  appears  to  be  poor,  for  two  reasons.  First,  not  all  the  lanes 
pass  close  to  the  seeded  items:  some  are  far  enough  to  measure  a  very  weak  field  and  yet  close 
enough  for  this  field  to  be  dominant  and  to  drive  the  inversion  algorithm.  In  these  cases,  the 
algorithm  tries  to  detect  the  target  despite  the  weak  measured  fields,  resulting  in  inaccurate 
estimates.  For  example  we  have  confirmed  that  the  first  shotput  (relative  Easting  about  2.75  m, 
relative  Northing  about  7  m)  is  properly  located  when  the  Meta  I  Map  per  travels  along  lanes  four 
and  five  from  the  left,  adjacent  to  the  target,  but  that  the  predictions  are  worst  along  lanes  three 
and  six.  The  associated  predictions  correspond  to  the  series  of  picks  between  3  m  and  4  m  in 
relative  Easting,  away  from  the  red  marker.  Flence,  lane  6  for  example  is  too  far  from  the  actual 
target  to  detect  it  accurately  and  given  the  noisy  data  and  weak  fields,  the  predictions  are  the 
best  the  algorithm  could  achieve.  The  second  reason  for  possible  mismatch  between  the  white 
crosses  and  the  field  peaks  in  Figure  4.2  is  that  the  picks  are  not  post- processed  in  the  figure 
and  not  averaged,  but  are  presented  as  raw  results.  We  again  refer  to  Figure  4.3  to  illustrate 
this  concept:  once  such  a  pattern  in  inverted  (x,y,z)  has  been  identified,  we  select  a  window 
of  detection,  where  the  target  is  believed  to  be  properly  picked.  In  this  case,  such  window  could 
be  defined  between  31  m  and  32  m.  Un-filtered  results  would  then  plot  all  the  picks  within  this 
window,  yielding  all  the  white  arrows  in  Figure  4.2.  A  post-processing  of  these  data  could  limit 
the  number  of  picks  to  the  point  where  y  =  0  (at  about  31.5  m  in  this  example),  or  average  all 
the  picks  within  31  m  and  32  m,  both  of  which  would  produce  a  single  white  pick.  Such  procedure 
is  adopted  in  the  next  example.  By  showing  unprocessed  picks  here,  however,  we  confirm  that 
the  algorithm  does  consistently  converge  to  reasonable  estimates  for  all  positions  that  are  close 
to  a  given  target,  thus  confirming  that  a  target  can  be  flagged  without  having  to  be  right  under 
the  sensor.  This  implies  that  the  MetalMapper  can  command  a  vehicle  to  stop  while  still  ahead 
of  a  target,  which  can  be  an  important  safety  feature  in  real  field  operations. 


The  second  example  pertains  to  a  similar  situation  with  data  collected  at  former  Camp  San 
Luis  Obispo  (SLO),  along  ten  lanes  spanning  an  area  of  about  7  m  by  60  m.  Unlike  Camp 
Butner,  the  subsurface  at  SLO  contains  much  less  clutter  items  as  can  be  seen  by  the  more 
uniform  background  in  Figure  4.4  (the  background  again  corresponds  to  the  interpolated  Hz 
value  at  the  central  receiver).  The  dynamic  data  were  post- processed  by  G&G  Sciences  Inc. 
using  a  method  based  on  matched  filter  to  a  sphere  response  [George  11],  thus  allowing  a  direct 
result  comparison  with  our  KF-EKF  approach. Results  are  shown  in  Figure  4.4  comparing  the  blue 
circles  (picks  by  G&G  Sciences  Inc.)  and  the  red  crosses  (averaged  picks  by  the  present  KF-EKF 
algorithm).  Like  in  the  previous  example,  our  target  picks  are  identified  by  regions  where  the 
inverted  (x,  y,  z)  positions  follow  trends  as  shown  in  Figure  4.3:  a  y  coordinate  matching  that  of 
the  moving  platform,  and  quasi-constant  x  and  z  coordinates.  Unlike  Figure  4.2,  however,  the 
range  of  values  corresponding  to  a  likely  target  are  here  averaged  to  produce  a  single  pick  on  the 
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Figure  4.2:  Target  picks  from  dynamic  data  collected  at  Camp  Butner,  NC.  The  background 
represents  the  interpolated  log10  values  of  the  Hz  component  at  the  central  receiver.  The  color 
scale  is  chosen  to  produce  red  and  blue  regions  at  location  of  strong  and  weak  magnetic  field 
values,  respectively. 
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Northing  [m] 


Figure  4.3:  Inverted  positions  around  a  Northing  position  of  31.5  m  in  Figure  4.2. 


map,  thus  resulting  in  much  less  red  cross  in  Figure  4.4  than  white  crosses  in  Figure  4.2.  The 
matching  with  results  from  G&G  Sciences  Inc.  is  seen  to  be  overall  very  good,  with  target  picks 
from  the  two  methods  overlapping  each  other.  Some  picks  from  G&G  Sciences  Inc.,  however, 
are  seen  to  stand  alone:  when  they  do  not  correspond  to  isolated  detection  points,  we  have 
confirmed  that  these  picks  could  be  matched  by  our  algorithm  by  relaxing  slightly  the  selection 
criterion.  For  example,  the  second  target  missed  on  lane  1  (dashed  left-most  lane)  could  have 
been  picked  if  we  had  allowed  variations  in  Easting  (x  coordinate)  and  depth  (z  coordinate)  of 
more  than  15  cm  within  the  detection  window. 

A  more  quantitative  comparison  between  the  inverted  positions  is  shown  in  Table  4.1 
where  we  also  compare  the  KF-EKF  predictions  with  those  of  our  reference  Gauss-Newton  algo¬ 
rithm  [Grzegorczyk  11],  Somewhat  expectedly,  it  is  seen  that  the  largest  error  corresponds  to 
the  component  in  the  traveling  direction  y  whereas  predictions  in  Easting  and  depth  are  typically 
within  a  few  centimeters.  The  larger  standard  deviations  are  due  to  the  fact  that  the  entire  de¬ 
tection  window  is  used  to  compile  these  results,  including  when  the  target  is  outside  the  physical 
range  of  the  sensor.  Although  we  have  confirmed  that  targets  do  not  need  to  to  be  exactly  under 
the  sensor  to  be  detected,  these  extreme  positions  typically  result  in  larger  errors  and  produce 
the  standard  deviations  between  10  cm  and  20  cm  reported  in  Table  4.1. 

As  a  further  validation,  we  compare  the  integrated  polarizabilities  obtained  by  G&G  Sciences 
Inc.  and  by  the  KF-EKF  algorithm.  Although  a  2.7  ms  data  collection  window  is  of  little  use  for 
classification  purposes,  it  can  provide  first  hand  volumetric  information  about  the  targets  that 
can  be  used  in  a  dig/no-dig  decision  making.  The  metric  is  the  root  mean  square  RMS^  of  the 
time  dependent  polarizabilities  across  the  N  time  gates  (between  times  t±  and  tN )  and  across  the 
entire  M- point  detection  window  (in  the  case  of  Figure  4.3  for  example,  the  detection  window 
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Mean  [cm]  Standard  deviation  [cm] 
\dx\  \dy\  \dz\  \dx\  \dy\  \dz\ 


KF-90ms  vs.  GN  0.5  4.4  1.3  9.6  18.2  6.9 

KF-90ms  vs.  G&G  0.7  14.1  0.6  9.2  15.7  8.7 


Table  4.1:  Mean  and  standard  deviations  (in  centimeters)  of  the  positions  obtained  by  (1) 
the  KF-EKF  algorithm  compared  to  the  Gauss-Newton  (GN)  reference  and  (2)  the  KF-EKF 
algorithm  compared  to  the  results  provided  by  G&G  Sciences  Inc..  The  sample  size  is  of  55 
targets  over  10  dynamic  lanes  at  Camp  San  Luis  Opbispo,  each  lane  being  about  60  m  in  length. 


is  defined  as  between  31  m  and  32  m  in  Northing  coordinates)  defined  as 


(4.1) 


Identifying  a  one  to  one  correspondence  between  G&G  Sciences  Inc.  and  the  KF-EKF  target 
picks  allows  a  one  to  one  comparison  between  the  corresponding  normalized  RMS^,  as  shown  in 
Figure  4.5.  Despite  magnitude  differences  of  the  three  larger  targets  (which  is  a  topic  of  further 
investigation  requiring  additional  data  collection  with  ground  truth),  all  the  most  prominent 
targets  are  identically  selected. 

3.  Note  on  processing  time 

Compared  to  the  previous  task,  our  purpose  was  to  speed  up  the  algorithm  and  achieve  “real 
time”  computational  efficiency.  This  was  motivated  by  a  series  of  discussions  while  working  with 
the  MetalMapper  in  a  field  setting,  where  the  need  of  a  “usable  algorithm”  was  emphasized,  i.e. 
a  code  that  can  actually  be  ported  onto  a  sensor's  control  panel  and  used  during  a  data  acquisi¬ 
tion  campaign.  We  therefore  reimplemented  our  algorithm  from  scratches,  leaving  the  KF  part 
(Kalman  filter)  essentially  unchanged  but  drastically  optimizing  the  EKF  part  (Extended  Kalman 
filter)  for  position  estimation.  As  a  result,  our  algorithm  currently  runs  in  about  250  ms  and 
provides  inversion  information  including,  location,  depth,  polarizabilities  (over  the  short  time  inter¬ 
vals  during  which  data  are  collected  in  a  dynamic  setting),  volume  information,  and  detectability 
criterion.  Various  metrics  have  been  tested  for  the  detectability  criterion  and  the  one  currently 
retained  is  based  on  a  linear  regression  of  the  inverted  polarizabilities  in  log-log  space. 

The  current  time  of  250  ms  is  obtained  on  a  regular  2  x  2.8  GHz  Quad-Core  Intel  Xeon 
computer.  Although  it  is  likely  to  be  a  more  powerful  computer  than  the  available  DAQ  currently 
associated  with  the  MetalMapper,  it  is  also  not  unrealistic  to  have  such  computer  power  on-board. 
In  addition,  we  have  not  performed  any  hardware  acceleration  such  as  using  GPU  or  parallelization. 
Hence,  if  necessary,  it  appears  that  it  would  be  very  possible  to  achieve  a  computation  time  of 
100  ms  as  we  had  originally  targeted,  and  realistically  call  our  algorithm  “real  time” . 

In  fact,  shorter  times  can  already  be  achieved:  the  KF-EKF  algorithm  runs  with  typically  ten 
realizations  and  loops  through  all  time  channels.  By  further  reducing  these  parameters,  we  can 
achieve  a  computational  time  of  about  70  ms  (hence  below  our  self-suggested  100  ms  limit),  yet 
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0  2  4  6 

Relative  Easting  [m] 

Figure  4.4:  Comparison  between  G&G  Sciences  Inc. target  picks  (blue  circles)  and  KF-EKF 
algorithm  target  picks  (red  crosses)  for  Camp  San  Luis  Obispo  dynamic  Meta  I  Map  per  data. 


at  the  expense  of  result  accuracy.  The  latter  is  being  currently  evaluated  and  preliminary  results 
are  shown  in  Figure  4.6.  The  figure  shows  a  history  of  estimated  positions  as  the  MetalMapper 
travels  by  a  target.  The  three  sub-panels  correspond  to  three  settings  of  the  algorithms,  from 
many  realizations  and  many  iterations  yielding  a  processing  time  of  850  ms,  to  fewer  realizations 
and  iterations  yielding  a  processing  time  as  low  as  70  ms.  In  all  cases  the  target  is  properly 
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Figure  4.5:  Comparison  between  the  averaged  polarizabilities  obtained  by  G&G  Sciences  Inc. 
(blue  line)  and  by  the  KF-EKF  algorithm  (red  line)  across  all  target  picks. 


located  and  locked  on,  but  with  varying  degrees  of  error.  A  quantitative  study  will  be  undertaken 
to  document  these  errors  more  specifically. 

Finally,  we  would  like  to  emphasize  that  the  KF-EKF  is  probably  best  applied  to  this  precise 
situation,  whereby  data  are  gathered  on  the  fly  over  targets  that  need  to  be  locked  on  for  a  certain 
time  period  (during  which  the  target  is  visible  to  the  Meta  I  Map  per).  We  put  this  situation  in 
contrast  to  waving  the  MPV  in  search  of  a  single  target:  in  this  case,  best  results  are  obtained 
when  all  data  points  are  used  simultaneously  in  order  to  obtain  good  inversion  accuracy.  This  is 
the  approach  we  adopt  in  our  Gauss-Newton  algorithm.  With  a  driving  MetalMapper,  however, 
the  problem  is  really  sequential:  each  shot  constitutes  a  sequential  information  that  can  be  used 
to  lock  onto  a  target  (as  it  passes  by  the  sensor),  or  unlock  from  a  target  that  moved  away  from 
the  sensor  field  of  view. 
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(a)  Setup  1:  850  ms  processing. 


(b)  Setup  2:  200  ms  processing. 
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Figure  4.6:  Newly  implemented  KF-EKF  algorithm  with  three  different  settings,  yielding  three  dif¬ 
ferent  data  processing  times.  The  figures  show  a  history  of  estimated  positions  as  the  MetalMap- 
per  drives  by  a  target. 
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(c)  Setup  3:  70  ms  processing. 
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4.  Graphical  User  Interface  proposal 

The  KF-EKF  algorithm  is  more  than  just  a  detection  algorithm  since  it  provides  an  estimate 
of  all  the  parameters  in  the  dipole  model  of  Eq.  (1.1):  position,  depth,  angles,  and  polarizabil¬ 
ities.  Moreover,  the  inverted  parameters  themselves  can  be  combined  to  extrapolate  additional 
information: 

(i)  Polarizabilities  can  serve  to  estimate  target  volume  and  three-dimensional  aspect  ratios, 

(ii)  Positions  and  depths  can  be  examined  for  consistency  at  sequential  sensor  positions  to 

inform  on  whether  the  target  is  likely  a  UXO  or  not, 

(iii)  Polarizabilities  can  be  matched  to  a  library  to  further  estimate  the  probability  of  UXO, 

(iv)  Further  processing  can  use  (F',0,0,  /3(f))  collaboratively  to  refine  all  estimates  mentioned 

above. 

The  information  content  is  therefore  rich  and  needs  to  be  displayed  appropriately  if  it  is  to  be 
used  by  an  operator  in  a  real-time  fashion.  An  example  of  a  graphical  user  interface  (GUI)  is 
shown  in  Figure  4.7  where  all  the  information  mentioned  above  is  displayed  in  real  time: 

>-  Left  column:  Position  information 

-  Top:  lateral  position  in  the  (xy)  plane  showing  the  inverted  values,  a  schematic  represen¬ 
tation  of  the  sensor,  and  the  field  background  (here  interpolated  between  all  receivers) 

-  Bottom:  Depth  as  function  of  past  sensor  positions,  which  allows  to  quickly  check  for 
consistency  (usually  associated  with  a  real  target  as  opposed  to  a  clutter  item) 

>-  Right  column:  Information  from  polarizabilities 

-  Top:  aspect  ratio  in  the  (xy)  and  (xz)  plane,  showing  an  ellipse  with  axes  (Px(t0),  /3y(t0)) 
and  (Px(to), Pz(,t0))  at  a  certain  time  t0. 

-  Middle:  Metric  indicating  the  probability  of  presence  of  a  target  of  interest.  In  this 
example,  the  metric  corresponds  to  a  correlation  of  the  three  polarizabilities  to  a  linear 
decay  within  the  first  1  ms  of  data  acquisition.  A  consistently  good  correlation  is  usually 
related  to  a  real  target  whereas  a  poor  correlation  is  associated  with  either  clutter  items 
or  algorithmic  divergences  (when  no  target  is  actually  present  beneath  the  sensor). 

-  Bottom:  Actual  inverted  polarizabilities  (/3x(t), Py(t), Pz(t))  here  shown  within  a  1  ms 
window. 

Figure  4.7  corresponds  to  a  snapshot  at  a  certain  position  of  the  MetalMapper  during  data 
acquisition  at  Camp  Butner.  The  full  movie  was  presented  along  with  a  summarizing  poster  at 
the  SERDP  Annual  Symposium  (Nov.  29  -  Dec.  1  2011,  Washington  DC). 
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KF-EKF  Real-time  Feedback  during  data  acquisition  campaign 


Real-time  detection  informatic 
as  the  sensor  travels  forwarc 
the  target  position  is  trackec 


^  Associated  depth  information 


Aspect  ratio  informs 
about  target  shape 


Likelihood  of  target 
detection  (threshold) 


Polarizabilities  (mostly  A 
volume  information)  ) 


Figure  4.7:  Snapshot  of  the  GUI  and  results  obtained  from  the  processing  of  dynamic  lane  33  at 
Camp  Butner.  Results  are  obtained  with  the  KF-EKF  algorithm  and  show  the  predicted  positions 
in  the  (xy)  plane,  the  predicted  depth,  the  inverted  polarizabilities  with  a  related  aspect  ratio  of 
the  UXO,  and  a  detection  metric  based  on  a  linear  fit  of  the  polarizabilities. 
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We  have  demonstrated  that  a  sequential  processing  of  EMI  dynamic  data  via  a  Kalman  filter 
(KF)  and  Extended  Kalman  filter  (EKF)  can  be  used  to  detect  UXO  and  give  good  estimates 
of  their  positions  and  polarizabilities.  The  demonstration  has  been  shown  with  data  from  two 
sensors,  the  MPV-II  and  the  MetalMapper,  in  very  different  settings  of  dynamic  data  collection. 
The  MPV-II,  meant  for  hand-held  operation,  was  used  to  interrogate  a  limited  area  within  which 
a  single  target  was  present,  akin  to  being  waved  by  an  operator  in  a  detection  mode.  The 
MetalMapper  instead  was  driven  along  adjacent  60-m  long  lanes  to  flag  target  locations  across 
a  large  geographical  area.  In  this  latter  configuration,  multiple  targets  were  present  in  the 
underground  and  sequentially  appeared  and  disappeared  from  the  field  of  view  of  the  sensor  as  it 
was  driven  across  the  area.  In  both  cases,  the  sensors  operated  in  dynamic  mode,  meant  for  rapid 
land  survey  at  the  expense  of  data  quality.  Despite  the  noisy  data,  the  KF-EKF  algorithm  was  able 
to  converge  to  proper  solutions  as  verified  independently  either  using  our  Gauss-Newton  reference 
algorithm,  using  ground  truth  information  whenever  available,  or  comparing  with  independently 
obtained  results. 

An  interesting  intrinsic  feature  of  the  KF-EKF  algorithm  is  its  sequential  processing  of  data 
as  they  become  available.  Plence,  at  every  new  data  collection  (e.g.  new  transmitter  position), 
the  algorithm  updates  the  estimated  parameters  based  on  the  newly  computed  Kalman  gain. 
Such  processing  can  be  performed  in  less  than  100  ms  on  a  regular  2  x  2.8  GHz  Quad-Core  Intel 
Xeon  computer,  which  is  the  real-time  limit  for  our  application.  The  inverted  parameters  include 
position,  depth,  and  polarizabilities,  albeit  limited  to  the  short  data  collection  time  range  inherent 
to  a  dynamic  survey.  Yet,  even  if  this  dynamic  data  is  not  sufficient  to  perform  classification,  it  is 
often  sufficient  to  identify  the  subsurface  anomaly  as  possible  UXO  or  not,  either  based  on  simple 
signal  correlation  to  a  library  or  based  on  a  volume  estimate.  The  method  can  therefore  be  used 
in  a  land  survey  for  real-time  target  mapping,  flagging  those  locations  to  which  the  instrument 
needs  to  return  for  more  exhaustive  cued  interrogation. 

A  current  limitation  of  our  algorithm  is  to  consider  the  process  of  Eq.  (2.1a)  as  a  Markov 
process  of  order  1,  i.e.  that  the  state  at  instance  k  depends  solely  on  the  state  at  instance  (k  —  1). 
For  UXO  detection,  however,  where  targets  are  by  definition  stationary,  this  can  be  a  limiting 
assumption.  For  example,  in  an  operating  mode  akin  to  that  of  the  MPV-II  in  Chapter  3,  multiple 
measurements  are  collected  atop  a  single  target.  The  polarizabilities  are  therefore  best  estimated 
by  using  all  data  points  simultaneously  in  the  inversion  algorithm.  This  is  the  approach  taken 
by  our  Gauss-Newton  algorithm,  whereby  data  from  all  grid  points  collectively  contribute  to  the 
construction  of  the  least-square  problem  from  which  a  single  position,  depth,  and  polarizability 
set  are  obtained.  The  assumption  of  a  Markov  process  of  order  1  is  therefore  a  restriction  to 
using  only  one  measurement  set  at  each  iteration,  whereas  all  previous  measurements  atop  the 
same  target  could  be  used.  In  a  setting  like  that  of  the  MetalMapper  in  Chapter  4,  however, 
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the  assumption  is  better  justified  since  the  sensor  drives  past  large  areas  and  is  not  continuously 
measuring  on  top  of  the  same  target.  We  therefore  recommend  the  use  of  the  present  KF-EKF 
algorithm  in  this  second  setting,  whereas  a  generalization  to  an  order  N  Markov  process  should 
be  implemented  for  the  inversion  of  data  collected  atop  a  single  target. 

In  addition,  several  avenues  of  follow  on  work  could  be  proposed  to  the  present  work: 

O  Integration  to  existing  systems:  The  present  KF-EKF  algorithm  is  currently  not  integrated 
with  any  sensor  but  stands  alone  and  performs  data  processing  off-line.  The  code  is  currently 
written  in  Matlab™  and  therefore  requires  the  acquisition  of  a  license  in  order  to  be 
executed.  In  order  to  perform  a  transfer  of  technology  and  incorporate  the  algorithm  directly 
onto  sensors  for  on-board  real-time  data  processing,  a  few  steps  need  to  be  undertaken: 
coordination  with  G&G  Sciences  Inc.  to  understand  platform  requirements,  transcription 
of  the  code  into  an  accepted  language,  on-board  testing  and  adjustments,  validation  with 
existing  data,  further  validation  during  data  acquisition. 

O  Real  field  validation:  Once  on-board,  the  algorithm  will  need  to  be  tested  on  a  per-sensor 
basis  (current  versions  have  been  validated  off-line  with  the  MPV-II  and  the  MetalMapper). 
Further  terrain  validation  with  these  sensors  will  be  necessary,  with  comparison  to  ground 
truth  whenever  possible,  or  comparison  with  different  methods  otherwise  (e.g.  our  Gauss- 
Newton  algorithm  or  others). 

O  Extension  to  new  sensors:  New  sensors  are  been  actively  developed  in  other  SERDP  pro¬ 
grams,  in  particular  the  Pedemis  [Barrowes  12]  which  is  meant  to  be  used  in  both  static 
and  dynamic  modes.  Owing  to  its  geometry,  the  Pedemis  can  be  viewed  as  an  extension  of 
the  MetalMapper  which  could  therefore  be  well  adapted  to  a  KF-EKF  real-time  processing 
(we  showed  in  this  work  that  a  dynamic  operation  of  the  MetalMapper  is  well  suited  with 
the  sequential  processing  feature  of  KF-EKF). 

O  Optimized  feedback  information:  The  GUI  proposed  in  Section  4.4.  is  merely  an  illustra¬ 
tion  of  a  possible  feedback  to  the  operator,  updated  in  real-time  as  data  are  processed 
sequentially.  The  proposed  GUI  not  only  displays  information  directly  obtained  from  the 
KF-EKF  algorithm  (position,  depth,  polarizabilities),  but  also  extrapolate  additional  infor¬ 
mation  such  as  volume  estimate,  aspect  ratio,  and  probability  of  being  a  target  of  interest 
based  on  a  certain  criterium.  The  latter  should  be  further  validated  with  real  data  in  order 
to  define  thresholds  of  probabilities  and  optimized  flagging  system  in  order  to  reduce  the 
rate  of  false  alarms.  The  GUI  should  also  be  written  in  a  language  accepted  by  the  sensors 
hardware. 

O  Longer  dynamic  collection  window:  The  KF-EKF  algorithm  is  more  advanced  than  a  sim¬ 
ple  detection  algorithm  since  it  provides  an  estimate  of  all  the  parameters  of  the  dipole 
model,  including  the  time  dependent  polarizabilities.  The  latter  are  currently  limited  to  a 
window  of  2.7  ms  which  is  too  short  to  attempt  a  reliable  target  identification.  As  a  result, 
the  polarizabilities  are  used  merely  to  estimate  target  volume  and  aspect  ratio.  With  little 
additional  complexity,  the  collection  window  could  be  extended  to  8  ms  and  the  polariz¬ 
ability  information  could  become  much  more  valuable,  opening  the  possibility  of  real-time 
target  identification  during  data  collection. 
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